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ABSTRACT 


Numerical  Simulation  of  Evaporating  Capillary  Jets 

Jason  D.  Zeda 
Advisor:  Dr.  Farzad  Mashayek 

A  detailed  numerical  study  of  evaporating  capillary  jets  is  presented.  The  analysis  is 
performed  through  use  of  a  Galerkin  finite  element  method  with  penalty  formulation  for  solving 
the  equations  of  motion  and  a  flux  method  for  tracking  the  free  surface. 

A  parametric  study  is  performed  to  analyze  the  temporal  instability  of  the  evaporating  jet. 
Through  varying  the  evaporation  rate,  Reynolds  number,  disturbance  wavenumber,  initial 
disturbance  amplitude,  and  density  ratio  the  outcomes  of  jet  breakup  are  investigated.  Also, 
pressure  distribution  inside  the  jet  and  multiple  satellite  drop  formations  are  analyzed.  Results 
are  compared  to  existing  analytical  conclusions  made  from  linear  stability  analysis. 

This  study  reveals  that  surface  evaporation  has  a  destabilizing  effect  for  the  low-speed 
jets,  which  are  considered  here.  That  is,  evaporation  flux  is  greater  at  the  neck  than  the  crest, 
which  accelerates  the  wave  growth.  Satellite  drops  also  reduce  in  size  as  evaporation  rate  is 
increased.  This  reduction  is  seen  in  both  the  radial  direction  due  to  vapor  leaving  the  surface  and 
along  the  axis  of  symmetry  due  to  decreased  breakup  time. 
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1.  Introduction 

The  instability  of  liquid  jets  has  been  a  classical  fluid  mechanics  problem  for  more  than  a 
century.  Early  investigations  focused  primarily  on  understanding  the  mechanisms  of  jet  breakup 
while  recent  jet  instability  studies  have  concentrated  on  more  accurate  predictions  of  growth  rate, 
drop  size  after  breakup,  and  effects  due  to  variations  in  initial  disturbance  amplitude,  disturbance 
type,  and  fluid  properties.  Despite  extensive  previous  studies  cited  in  the  literature  and  briefly 
reviewed  below,  it  has  been  only  recently  that  the  dynamics  of  a  jet  undergoing  surface 
evaporation  has  been  investigated  by  Lian  and  Reitz  (1993)  using  a  linear  stability  analysis. 
While  this  study  reveals  several  new  and  interesting  phenomena,  it  is,  understandably,  unable  to 
predict  accurate  results  near  the  breakup  point  due  to  the  assumption  of  infinitesimal  perturbation 
of  the  liquid  jet. 

The  objective  of  the  present  work  is  to  relax  some  of  the  assumptions  adopted  by  Lian 
and  Reitz  (1993)  and  to  provide  a  more  versatile  computational  analysis  of  evaporating  jets  by 
implementing  a  Galerkin  finite  element  method  in  conjunction  with  a  height  flux  method  (HFM) 
for  interface  tracking,  which  was  developed  by  Mashayek  and  Ashgriz  (1993).  The 
dimensionless  continuity  and  momentum  equations  are  solved  on  a  moving  mesh  consisting  of 
four-node  quadrilateral  elements.  Considering  an  axisymmetric  incompressible  Newtonian 
liquid  jet  of  infinite  length  initially  subjected  to  a  periodic  surface  disturbance  a  parametric  study 
is  performed  to  show  how  evaporation  affects  capillary  jet  instability  while  varying  liquid 
Reynolds  number  (Re),  disturbance  wavenumber  (k),  disturbance  amplitude  (£0),  evaporation  rate 
(P),  and  density  ratio  (s).  Similar  to  Lian  and  Reitz  we  focus  on  the  temporal  evolution  of  the  jet 
thus  only  one  wavelength  is  examined  at  a  time  -  in  practice,  we  solve  for  only  one  half  of  this 
wavelength  due  to  the  symmetry  of  the  jet.    Also  similar  to  Lian  and  Reitz  (1993),  the  surface 


evaporation  of  the  perturbed  jet  is  approximated  by  that  of  a  spherical  drop  with  the  identical 
curvature  as  the  local  surface  curvature  of  the  jet.  In  this  manner,  there  is  no  need  to  solve  for 
the  gas  phase. 

The  observation  that  a  liquid  jet  issuing  from  a  nozzle  will  eventually  breakup  into  small 
drops  when  subjected  to  even  minute  disturbances  led  Rayleigh  (1879)  to  provide  the  first 
analytical  description  for  the  temporal  instability  of  an  inviscid,  incompressible  jet  using  a 
normal  mode  analysis.  He  showed  that  an  axisymmetric  harmonic  disturbance  of  the  form: 

r\  =  e0exp(pt-ikz),  (1) 

would  grow  in  time  according  to: 


(0  = 


l^k){i-e)k112 


(2) 


where  co  is  the  growth  rate,  /  is  the  time,  and  /0  and  I\  are  the  modified  Bessel  functions  of  the 
first  kind.  In  (1)  and  (2),  length  and  time  variables  have  been  non-dimensionalized  by  r0  and 
fo/ve,  where  ve  is  the  capillary- wave  velocity  (a/pr0)'/2.  Here,  r0  is  the  undisturbed  radius  of  the 
jet,  a  is  the  liquid  surface  tension  coefficient,  and  p  is  the  density  of  the  liquid.  This  result 
predicts  that  disturbances  grow  in  time  for  k<\  and  oscillate  for  k>\,  and  the  maximum  growth 
rate  occurs  at  £=0.697.  In  other  words,  if  the  wavelength  (k=2nrjk)  of  an  axisymmetrical 
disturbance  is  greater  than  the  cylinder  circumference  the  jet  will  be  unstable  and  if  the 
disturbance  wavelength  is  less  than  the  cylinder  circumference  the  jet  will  be  stable. 

For  over  a  century  since  Rayleigh  noticed  that  surface  tension  had  to  work  against  inertia 
in  causing  instability  researches  have  continued  to  extend  his  work  in  both  theoretical  and 
experimental  realms.  For  example,  viscous  effects  have  been  considered  by  Weber  (1931),  who 
analyzed  the  motion  of  thin  sections  of  fluids,  and  Chandrasekhar  (1961),  whose  book  focused 


solely  on  the  subject  of  instability  using  the  Navier-Stokes  equation  to  examine  it.  To  test  the 
results  of  linear  theory  quantitative  experiments  were  performed  by  such  researchers  as  Goedde 
and  Yuen  (1970).  The  experiments  clearly  identified  that  drop  size  varied  at  breakup  leading  to 
the  conclusion  that  nonlinear  behavior  became  dominant  during  instability.  Yuen  (1968)  was  the 
first  to  extend  Rayleigh's  linear  analysis  by  considering  a  third-order  perturbation  expansion, 
which  would  give  some  explanation  of  the  production  of  satellite  drops  between  the  larger  main 
drops.  Further  experimental  investigation  of  both  satellite  and  main  drops  have  been  done  by 
Vasallo  and  Ashgriz  (1991)  and  Kowalewski  (1996),  which  looked  at  the  production  of  multiple 
satellites  and  the  description  of  the  jet  surface  near  breakup,  respectively. 

Due  to  recent  development  of  computational  facilities  and  numerical  techniques  flow 
simulations  have  become  common  place.  Inviscid  flow  was  studied  by  Mansour  and  Lundgren 
(1990)  and  Stokes  flow  has  been  considered  by  Stone  and  Leal  (1989).  These  two  simplifying 
factors  breakdown  close  to  breakup  because  both  viscous  and  inertial  effects  become  important. 
In  this  case  the  full  Navier-Stokes  equation  has  to  be  considered  in  the  fluid  domain.  A  finite 
difference  one-dimensional  model  was  considered  by  Eggers  and  Dupont  (1994)  to  simulate  drop 
formation.  The  most  extensive  study  of  jet  breakup  using  the  Navier-Stokes  equation  was 
performed  by  Ashgriz  and  Mashayek  (1995)  using  a  Galerkin  finite  element  method.  For  a 
review  of  the  computational  methods  used  for  free  surface  flows  see  Tsai  and  Yue  (1996). 
Further  reviews  of  the  problem  of  jet  instability  are  given  by  Bogy  (1979)  and  more  recently  by 
Eggers  (1997). 

Unfortunately,  little  mention  is  made  about  evaporating  jets  when  studying  the  problem 
as  a  whole.  This  is  primarily  due  to  the  complexity  of  determining  an  evaporation  equation  for 


an  irregular  shaped  surface.    Hopefully  future  experimentation  will  be  performed  to  verify  the 
spherical  model  used  herein. 


2.  Formulation  and  Methodology 

Consider  a  viscous  liquid  circular  jet  injected  with  velocity  U  into  an  inviscid  gas  at 
temperature  7^,  in  zero  gravity.  Since  we  are  dealing  with  interactions  of  the  gas  with  only  one 
jet,  it  can  be  assumed  that  the  total  heat  capacity  of  the  gas  is  much  larger  than  that  of  the  liquid 
and  T^,  far  from  the  surface  of  the  jet,  is  constant.  The  implementation  of  an  evaporation  model, 
described  below,  eliminates  the  need  for  the  solution  of  the  energy  equation;  thus  the  governing 
equations  include  only  mass  and  momentum  conservations.  Assuming  both  the  liquid  and  the 
gas  are  incompressible  with  constant  properties,  these  equations  are  described  as: 

Vw,=0,  (3) 

^i  +  w    Vw,.=--  V^.+v.-W,    i  =  l,g  (4) 

dt  p, 

where  /  and  g  stand  for  liquid  and  gas,  respectively.  In  (3)  and  (4),  w,  p,  p,  and  v  denote  the 
velocity,  pressure,  density,  and  kinematic  viscosity,  respectively.  For  a  liquid  jet  evaporating 
with  a  flux  m ,  the  jump  condition  and  the  normal  stress  balance  on  the  interface  are  given  by 
(Lian  and  Reitz,  1993): 

th  =  P/  (w,  -  ws  )•  n  =  p^  (wg  -  w t  )■  n,  (5) 

m(wg  -  w,  )■  n  +  |n  •  (xg  -  x,  )J- n  +  oV  •  n  =  0,  (6) 

where  a,  x,  ws,  and  n  denote  the  surface  tension  coefficient,  the  stress  tensor  including  the 
pressure  component,  the  interface  velocity,  and  the  outward  unit  normal,  respectively.   The  first 


term  on  the  left-hand  side  of  (6)  is  the  recoil  force  on  the  jet  due  to  surface  evaporation.   This 
force  stems  from  a  net  momentum  that  is  impacting  the  surface.  Eliminating  ws  in  (5)  yields: 


(w,-w,)-n  = 

r 


p/ 


m 

1-  — 

I       ^  j 
In  this  study  we  consider  pg«P/,  therefore  the  effects  of  the  fluctuations  in  the  ambient  gas 

pressure  may  be  neglected  in  comparison  to  the  effects  of  the  pressure  fluctuations  in  the  liquid. 
Further,  the  mean  pressure  of  the  gas  affects  only  the  mean  pressure  of  the  jet  and,  for 
incompressible  liquid,  does  not  influence  the  instability  of  the  jet.  As  a  result,  we  can  neglect  the 
ambient  gas  pressure  and,  for  inviscid  gas,  eliminate  %g  from  (6).  Then,  by  implementing  (7),  (6) 
yields: 

n  x7   n  =  oV  n  +  — .  (8) 

P* 

This  simplifies  the  problem  significantly,  as  there  will  be  no  need  to  solve  for  the  gas  phase.  It 
must  be  noted  that,  in  this  manner  the  mean  pressure  of  the  jet  is  calculated  relative  to  the  mean 
pressure  of  the  ambient  gas. 

There  is  no  analytical  solution  for  evaporation  of  a  jet  with  a  perturbed  surface. 
However,  following  Lian  and  Reitz  (1993)  one  may  approximate  the  unit  area  evaporation  rate 
of  a  wavy  curved  surface  of  local  radius  of  curvature  Rc  with  that  of  a  sphere  of  the  same  radius. 
In  this  study,  assuming  that  the  surface  temperature  of  the  jet  is  always  at  the  boiling  temperature 
Tb,  the  flux  of  evaporation  is  modeled  as: 

m  =  Q-     for  R>0,  (9) 

Rc 

where, 


<H 


(  K     \ 

^|ln 
V     pg  J 


uCjT„-Tj 


(l  +  0.276  Re  /Vr/3),  (10) 


where  k,  Cp,  and  Lv  denote  the  thermal  conductivity,  the  specific  heat,  and  the  latent  heat  of 
evaporation,  respectively.  The  Reynolds  number,  Re,  and  the  Prandtl  number,  Pr,  are  defined  as 
Reg=2UrJvg  and  Prg=pgvgCpg/K.g.  Equations  (9)  and  (10)  are  a  modified  form  of  the  d2-law  for 
evaporation  of  a  single  drop  in  zero  gravity  (Williams  1985  and  Turns  1996).  The  modification, 
here,  is  by  substituting  the  local  radius  of  curvature  of  the  perturbed  jet  for  the  radius  of  the 
spherical  drop.  This  modification  was  first  introduced  by  Lian  and  Reitz  (1993).  The  physical 
reasoning  is  that  the  surface  of  the  deformed  jet  may  be  "locally"  considered  as  the  surface  of  a 
spherical  drop  having  the  same  radius  of  curvature.  It  is  noted  that  the  utilization  of  (9)  is 
limited  to  surface  deformations  such  that  Rc  >  0. 

In  the  remaining  equations  (for  the  liquid  phase)  an  axisymmetric  incompressible 
Newtonian  liquid  jet  in  zero  gravity  is  considered.  The  variables  are  non-dimensionalized  by  the 
undisturbed  jet  radius  r0,  and  the  characteristic  time  (p/"03/o)/2.  The  dimensionless  continuity  and 
momentum  equations  are: 

V-v  =  0,  (11) 

Re  — +  v-Vv  |=V-T,  (12) 

I* 

where  v=W/=(w,v)  is  the  velocity  vector  and  T=i/=-/?I+[Vv+(Vv)T]  is  the  stress  tensor  for  a 
Newtonian  fluid  with  p  now  denoting  the  normalized  liquid  pressure.  The  Reynolds  number  in 
(12)  is  defined  as  Re={\Ni){<5rJpi)v\  Non-dimensionalizing  (8)  the  stress  balance  on  the 
interface  yields: 

Tn  =  Re(K  +  aK2)n,  (13) 


where  K  is  the  curvature  of  the  surface  and 


S  = 


(14) 


As  a  result  of  the  non-dimesionalization,  the  parameters  affecting  the  instability  of  an 
evaporating  jet  reduce  to  Reynolds  number  Re,  the  dimensionless  evaporation  rate  p\  and  the 
density  ratio  s.  In  the  next  section  we  investigate  the  effects  of  these  parameters  as  well  as  the 
effects  of  initial  disturbance  amplitude  e0  and  wavenumber  k  via  numerical  simulation. 

In  this  investigation  it  is  assumed  that  the  free  surface  can  be  represented  by  a  height 
function  h(z,t),  as  shown  in  Fig.  1,  thus  the  curvature  K  is  calculated  using: 

h„  1 


K  = 


(15) 


(i  +  /cf    h{\  +  hlf 

Since  a  temporal  analysis  is  considered  symmetry  boundary  conditions  can  be  applied  on  the  axis 
and  planes  of  symmetry: 


r,  v 


Figure  1 :  Fluid  subvolumes  represented  by  their  height  and  thickness. 


du 


0   and    v  =  0,   at  r  =  0 


w  =  0   and   —  =  0,   at  z  =  0,  — , 

3z  2 


(16) 
(17) 


A  Galerkin  finite  element  method  is  used  to  solve  (11)  and  (12).  Using  a  penalty  function 
formulation  (Hughes,  Liu  and  Brooks  1979),  the  pressure  term  is  eliminated  from  the  set  of 
unknown  variables  by  absorbing  the  continuity  equation  into  the  momentum  equation.  In  this 
formulation  the  pressure  is  defined  as: 

P-YVv,  (18) 

where  Y  is  a  large  number  of  order  (109)  depending  on  viscosity  and  Reynolds  number.  Four- 
node  bilinear  isoparametric  elements  are  used  to  approximate  the  velocity  distribution  over  each 
element: 


\(r,z,t)  =  ^\,(t)Ni(r,z,t). 


(19) 


i=i 


A  moving  mesh  is  considered  to  discretize  the  computational  domain  making  the  shape 
functions,  Nj,  time  dependent.  In  order  to  obtain  the  finite  element  formulation,  (12)  is 
multiplied  by  the  shape  function  Nj  and  integration  is  carried  over  the  element  volume.  After  the 
divergence  theorem  is  invoked  the  resulting  equation  is: 


J  k 

Jam  I     J 


Re 


9v 
dt 


+  v-Vv 


A 


J 


+  VN]  -[-/?I  +  |yv  +  (Vv)r]  \dQ=  \    NjJ-ndT, 


(20) 


where  Q  is  the  volume  and  T  is  the  surface  area  of  the  element.  Through  substitution  of  (13)  and 
(18)  into  (20)  we  obtain  the  following  closed  form  finite  element  formulation: 


f     \NjRe 


dt 


+  v- Vv 


+  VN]    ^(Vv^  +  [vv  +  (Vv)7]|JQ  =  J    N^K  +  aK^dT.  (21) 


The  above  formulation  is  based  on  a  fixed  mesh.  In  the  present  problem  the  mesh  is  moving, 
therefore  special  treatment  of  the  time  derivatives  is  required.  Since  the  shape  function  is  time 
dependent,  the  time  derivative  of  velocity  in  discretized  form  becomes: 


=  /  — LNi  +  >  v,.  — - 


(22) 


V  w    Jr.Z  '=1      **»  '=1  V  JTA 

The  last  term  of  (22)  introduces  a  new  convective  term  in  (21).  Here,  we  allow  the  motion  of  the 

nodes  in  the  r-direction  only,  according  to  the  following  simple  rule: 

z;(?+8f)=z,(0=constant,    r,(f+8r)=cr((0,  (23) 

where  the  subscript  i  refers  to  the  node  number,  and  c=c(z,t)  is  a  constant  for  each  column  of 

nodes  in  the  radial  direction  defined  as: 

c  =  h(z,t  +  5t)/h(z,t).  (24) 

In  this  case,  Mashayek  and  Ashgriz  (1993)  have  shown  that  the  total  derivative  of  velocity 

becomes: 

D\  _  3v     (      c -1    }d\        d\ 
Dt  ~  dt 


v r  | —  +  u — .  (25) 

^       bt     jdr        dz 


Substitution  of  (25)  into  (21)  completes  the  finite-element  formulation  for  the  present  problem. 

The  free  surface  of  the  jet  is  determined  using  the  Height  Flux  Method  (HFM)  developed 
by  Mashayek  and  Ashgriz  (1993).  As  shown  in  Fig.  1  the  domain  is  divided  into  several  vertical 
subvolumes  of  width  6zi  and  volume  V/.  The  location  of  the  free  surface  on  the  left  and  right 
hand  sides  of  a  subvolume  is  given  by  /z,  and  hi+\  at  time  t,  respectively.  At  the  end  of  each  time 
step  fluxes  of  the  fluid  from  each  subvolume  to  its  neighboring  subvolumes  are  calculated  using 
the  velocity  field  determined  by  the  finite  element  solution  of  the  governing  equations.  Next, 
with  the  knowledge  of  the  evaporation  rate,  the  volume  change  bVe  of  each  subvolume  due  to 
evaporation  at  the  interface  is  calculated  using: 


8V  =  2ti 


(h^+hA 


M- '-  1^(35/5/, 

2 


(26) 


where  5/  indicates  the  length  of  the  interface  confined  within  the  subvolume.  The  volume  of 
fluid  within  each  subvolume  is  further  modified  by  calculating  the  flow  through  the  side  planes. 
Then,  it  is  assumed  that  the  part  of  the  interface  located  between  any  neighboring  pair  of 
subvolumes  can  be  approximated  by  a  line  segment  described  as  h=az+b,  where  a  and  b  are 
constants  which  can  be  determined  by  the  known  volumes.  More  details  and  the  accuracy  of  this 
method  are  discussed  by  Mashayek  and  Ashgriz  (1993). 


Mesh  Size 

4x30 

4x40 

4x50 

4x60 

Growth  Rate 

0.712 

0.712 

0.712 

0.712 

Table  1:  Typical  convergence  test  for  growth  rate:  7te=200,  eo=0.05,  k=0.1,  P=3xl0"2,  and  s=\0 
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Before  numerical  results  were  compiled  a  convergence  study  was  performed  to  determine 
the  effect  of  mesh  size  on  growth  rate.  The  results  are  shown  in  table  1.  From  this  we  conclude 
that  grid  size  does  not  affect  growth  rate  since  only  the  height  of  the  free  surface  at  the  neck  and 
swell  where  of  importance  for  calculating  this  value.  To  save  computational  time  the  majority  of 
cases  were  ran  with  a  minimal  number  of  elements.  The  mesh  size  for  the  determination  of 
growth  rate  was  limited  to  four  elements  in  the  radial  direction  and  the  axial  elements  were 
changed  based  on  the  wavenumber  as  given  in  table  2.  For  the  cases  where  pressure  distribution 
and  velocity  vectors  were  plotted  the  grid  was  refined  in  order  to  obtain  better  print  quality  and 
color  dispersion. 


Wavenumber 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

Number  of  elements 

50 

40 

40 

30 

30 

20 

20 

Table  2:  Number  of  axial  elements  used  in  simulations  for  determination  of  growth  rate. 
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3.  Effect  of  Evaporation  on  Jet  Instability 

We  consider  the  temporal  evolution  of  a  liquid  jet,  initially  at  rest  with  a  spatially 
harmonic  surface  described  as: 

r=R-e0cos(kz),  (27) 


where  R  is  determined  such  that  the  volume  of  the  jet  is  kept  constant  when  the  initial  amplitude 
is  changed: 

R  =  ^-i<J2-  (28) 

Due  to  symmetry,  only  half  a  wavelength  is  considered.  The  trough  (neck)  of  the  initial  surface 
is  set  at  z=0  and  the  crest  (swell)  at  z=A72.  The  dynamics  of  this  jet  due  to  capillary  and  recoil 
forces  are  investigated  for  various  values  of  dimensionless  evaporation  rate,  density  ratio, 
Reynolds  number,  disturbance  wavenumber,  and  initial  disturbance  amplitude. 

Detailed  description  of  the  shape  evolution  of  liquid  jets  with  Re=200  and  10  is  presented 
in  Figs.  2  and  3  for  £=0.7,0.4  and  (3=0,  10"2,  and  3xl0"2.  The  following  characteristics  for  the 
breakup  of  an  evaporating  capillary  jet  can  be  observed  from  Figs.  2  and  3:  (i)  As  evaporation 
rate  increases  the  breakup  time  decreases;  (ii)  As  Reynolds  number  decreases  the  breakup  time 
increases;  (iii)  The  length  and  diameter  of  the  liquid  ligament  decreases  with  increasing 
evaporation  rate  and  decreasing  Reynolds  number;  (iv)  The  breakup  point  is  closer  to  the  initial 
trough  as  evaporation  rate  increases,  due  to  the  decrease  in  time  that  nonlinear  effects  can 
influence  jet  breakup.  A  more  detailed  description  of  these  characteristics  is  given  in  the 
following  sections. 
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Figure  2:  Time  evolution  of  the  instability,  /te=200,  eo=0.5,  s=10~3:  (a)  k=0J,  (3=0;  (b)  k=0J, 
P=102;  (c)  k=0.1,  p=3xl0"2;  (d)  k=0A,  p=0;  (e)  k=0A,  (3=102;  if)  Jfc=0.4,  (3=3xlO'2. 


3.1  Location  of  the  breakup  point 

Linear  theory  predicts  that  the  breakup  point  is  always  at  the  trough  of  the  initial 
disturbance  where  the  initial  jet  radius  is  minimum.  However,  Figs.  2  and  3  clearly  show  that  the 
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Figure  3:  Time  evolution  of  the  instability,  Re=\0,  eo=0.5,  s=\0'3:  (a)  k=0.1,  (3=0;  (b)  fc=0.7, 
(3=10~2;  (c)  £=0.7,  (3=3xl02;  (d)  £=0.4,  (3=0;  (e)  £=0.4,  (3=10~2;  (/)  £=0.4,  (3=3xl02. 


breakup  point  moves  in  the  direction  of  the  swell  towards  the  end  of  the  simulation.  The 
capillary  forces  on  the  jet  tend  to  "squeeze"  the  liquid  in  the  neck  region  into  the  swell.  This 
results  in  the  formation  of  two  distinct  regions  in  the  jet.  The  middle  region  where  the  curvature 
in  the  r-z  plane  is  very  small,  resembling  a  liquid  ligament,  and  the  swell  region  where  the 
curvature  in  the  r-z  plane  is  large.  Note  from  Figs.  2  and  3  that  the  disturbances  stay  sinusoidal 
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for  a  very  long  time  (relative  to  the  breakup  time)  and  agree  with  linear  theory  until  non-linear 
effects  begin  to  move  the  location  of  the  minimum  radius. 

In  order  to  gain  a  better  understanding  of  the  development  of  the  breakup  point  of  the  jet, 
we  have  plotted  the  motion  of  the  minimum  radius  along  the  jet  from  time  t=0  to  the  breakup 
time.  This  plot  is  shown  in  Fig.  4  for  a  jet  with  Re=200,  eo=0.05,  s=10"3,  and  for  three  different 
evaporation  rates.  Two  distinct  time  periods  can  be  identified  for  each  of  the  curves  in  these 
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Figure  4:  Motion  of  the  minimum  radius  along  the  jet  axis:  (a)  (3=10  ;  (b)  P=2xl0  ;  (c)  P=3xl0 


14 


figures.  The  most  important  being  the  time  that  the  minimum  point  stays  at  z=0.  This  time 
represents  over  90%  of  the  total  breakup  time  in  all  cases  and  it  first  reduces  and  then  increases 
with  an  increase  in  the  wavenumber.  The  case  with  the  maximum  disturbance  growth  rate  is  the 
first  to  begin  movement  away  from  z=0.  The  second  distinct  time  period  represents  the  rapid 
movement  of  the  minimum  jet  radius  from  the  neck  towards  the  swell  region.  This  sudden 
relocation  of  minimum  point  in  the  axial  direction  is  larger  for  smaller  wavenumbers.  The  end 
points  on  all  the  curves  in  Fig.  4  represent  the  breakup  times.  There  are  two  interesting  features 
to  note  from  Fig.  4.  The  first  is  that  as  evaporation  rate  increases  the  location  of  the  minimum 
point  is  closer  to  the  neck  for  each  wavenumber.  This  represents  the  production  of  a  smaller 
middle  ligament  or  satellite  drop.  The  second  feature  reveals  that  as  evaporation  rate  increases 
the  range  of  time  that  each  wavenumber  breaks  away  from  z=0  reduces. 


3.2  Growth  rate  of  the  disturbances 

Linear  theory  not  only  defines  the  region  of  unstable  disturbance  wavenumbers  but  also 
provides  their  growth  rates.  These  growth  rates  are  useful  in  estimating  the  breakup  time  and 
length.  According  to  linear  theory  the  variation  of  the  logarithmic  value  of  the  amplitude  of  the 
surface  disturbances  with  time  is  linear.  However,  for  an  actual  liquid  jet  this  amplitude 
variation  is  nonlinear  at  the  beginning  of  the  simulation  due  to  the  fact  that  radial  velocity  of  the 
jet  starts  from  zero.  The  amplitude  variation  is  nonlinear  also  near  the  breakup  point  due  to  the 
movement  of  the  minimum  radius.  For  a  selected  few  of  our  simulations  we  have  plotted  the 
logarithmic  values  of  the  normalized  amplitudes  of  the  neck  ((/?e-rn)/£o),  swell  ((rs-Re)/E0),  and 
their  difference  ((rs-rn)/e0)  in  Fig.  5.  Here,  rn  and  rs  are  the  radii  of  the  neck  and  swell, 
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Figure  5:  Variation  of  the  amplitude  of  the  neck  (black),  swell  (red),  and  their  logarithmic 
difference  (green),  Re=200,  £o=0.5,  s=10"3:  (a)  k=0A,  p=10"2;  (b)  £=0.7,  p=10  2;  (c)  A:=0.4, 
P=2xl0"2;  {d)  £=0.7,  p=2xl0~2;  (e)  £=0.4,  p=3xl  0"2;  (/)  £=0.7,  [3=3x1 0"2. 


respectively,  and  Re  is  the  instantaneous  radius  of  the  equivalent  cylindrical  jet,  which  changes 
during  the  simulation  due  to  evaporation.    Notice  in  Fig.  5  that  the  logarithmic  variation  of 
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amplitude  of  the  neck  or  swell  is  not  perfectly  linear.  Figure  5  reveals  interesting  features  for  the 
temporal  evolution  of  amplitudes  of  the  neck  and  the  swell  with  respect  to  each  other.  Initially 
the  growth  rate  of  the  neck  point  is  larger  than  the  growth  rate  of  the  swell  point,  since  a  radial 
displacement  in  the  neck  region  corresponds  to  a  smaller  radial  displacement  in  the  swell  region 
for  the  same  volume  displacement.  Due  to  the  formation  of  a  liquid  ligament,  i.e.  minimum 
radius  point  moving  away  from  the  neck,  the  growth  rate  of  the  amplitude  of  the  neck  decreases 
eventually  dropping  below  that  of  the  swell  for  the  longer  wavelengths.  This  is  clearly  shown  in 
Fig.  5a,c,e,  which  also  indicates  that  as  we  increase  the  evaporation  rate  the  time  after  the  cross 
point  is  reduced.  This  represents  the  formation  of  a  smaller  liquid  ligament,  which  can  be  seen 
in  Figs.  2  and  3.  Figure  5b,d,f  is  representative  of  smaller  disturbance  wavelengths.  Due  to  the 
reduced  time  of  breakup  non-linear  effects  are  only  visible  in  the  last  instances  of  instability. 
Comparing  Fig.  5a,c,e  with  5b,d,f  we  see  that  the  shape  of  the  difference  curve  varies  only  in 
slope  when  a  change  in  wavenumber  is  made. 

Similar  to  previous  research  in  this  area  we  use  the  linear  region  of  the  difference  curve 
to  provide  values  that  can  represent  the  growth  rate  for  each  wavenumber.  Therefore,  a  line  is  fit 
to  the  difference  between  the  amplitudes  of  the  neck  and  swell  and  its  slope  is  measured.  The 
line  is  fit  only  to  the  middle  region  of  each  curve  and  the  end  portions  are  ignored.  The 
dispersion  curves  are  obtained  using  the  calculated  values  of  the  growth  rate  for  different 
wavenumbers  and  for  various  values  of  p\  Re,  e0,  and  s.  The  data  are  plotted  in  Fig.  6.  As 
mentioned  earlier  Fig.  6a  shows  quantitatively  that  as  we  increase  the  evaporation  rate  the 
growth  rate  increases  therefore  reducing  the  time  it  takes  for  the  jet  to  breakup  into  droplets.  All 
curves  in  Fig.  6a  give  the  maximum  growth  rate  at  k=0.7,  similar  to  linear  theory.  Figure  6b 
presents  the  effects  of  Reynolds  number  on  growth  rate  of  disturbances.  Here  we  can  clearly 
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Figure  6:  Growth  rate  of  disturbances:  {a)  Re=200,eo=0.05,s=l0~3;  (b)  eo=0.05,p=3xl0"2^=10"3; 
(c)  /te=200,p=3x!0  2,s=W3;  (d)  /?e=200,eo=0.05,p=3xl0"2;  (e)  Lian  &  Reitz  (1993). 


identify  for  evaporating  jets  that  as  we  decrease  the  Reynolds  number  the  growth  rate  of 
disturbances  decreases  due  to  the  more  viscous  nature  of  the  jet  not  allowing  the  fluid  to  be 
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displaced  as  easily.  Also  from  this  graph  we  see  a  movement  of  the  maximum  growth  to  a  lower 
wavenumber  for  the  case  of  highly  viscous  fluids.  Figure  6c  demonstrates  the  effects  of  initial 
disturbance  amplitude  on  growth  rates.  As  mentioned  earlier,  linear  theory  assumes  infinitesimal 
initial  disturbance  amplitude.  Therefore,  for  large  initial  disturbance  amplitudes  the  theoretical 
solutions  result  in  large  errors.  We  see  from  Fig.  6c  that  as  we  increase  this  value  growth  rate 
decreases  at  each  wavenumber.  However,  this  does  not  mean  that  the  jet  will  take  longer  to 
breakup  for  large  initial  disturbance  amplitudes.  This  figure  only  gives  meaning  to  the  rate  and 
not  the  overall  distance  that  the  radius  has  to  decrease  in  order  to  pinch  off  the  jet.  Figure  6d 
provides  insight  into  the  effect  of  the  density  ratio  on  growth  rate.  Here  we  clearly  see  that  as  we 
increase  the  density  of  the  gas  or  decrease  the  density  of  the  liquid  (i.e.  larger  s)  the  growth  rate 
of  the  evaporating  jet  decreases.  This  can  be  explained  through  the  reduction  of  the  recoil  force, 
where  the  net  momentum  imparted  on  the  surface  is  diminished  by  increasing  s  (see  (13)). 

Figure  6e  shows  the  results  from  Lian  and  Reitz's  (1993)  linear  stability  analysis  for  an 
evaporating  viscous  jet  with  s=0.25xl0"8,  We^lO"4,  and  Z=1CT4,  where  We^  represents  the 
Weber  number  and  Z  represents  the  Ohnesorge  number.  Since  we  do  not  solve  for  the  gas  phase 
we  must  convert  their  formulation  to  agree  with  our  input  parameters.  For  the  case  shown  in  Fig. 
6e  this  would  result  in  a  /te=20,000  due  mainly  to  their  very  small  value  for  density  ratio  (e.g.  s 
for  air  and  water  would  be  approximately  1.2x10"  ).  Using  this  Reynolds  number  and  similar 
values  of  s  and  (3  we  attempted  to  correlate  a  direct  comparison  but  was  unable  to  reproduce  their 
results.  The  calculated  growth  rates  in  our  simulations  were  very  similar  to  cases  with  no 
evaporation.  Only  when  we  increased  (3  to  around  5xl0"5  did  we  obtain  growth  rates  larger  than 
that  of  non-evaporating  cases.  Although  we  could  not  reproduce  the  exact  results  from  the 
analytical  study,  even  using  the  same  parameters,  we  both  conclude  the  same  general  trends  for 
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growth  rate.  As  evaporation  rate  increases  growth  rate  increases  for  every  wavenumber.  One 
possible  reason  for  the  discrepancies  between  our  results  and  theirs  is  the  linear  theory's 
assumption  of  infinitesimal  disturbance  amplitude.  For  non-evaporating  jets  the  growth  rate 
does  not  vary  up  to  an  initial  disturbance  amplitude  of  approximately  0.05.  However,  Fig.  6c 
clearly  shows  that  the  growth  rate  varies  in  this  region  for  evaporating  jets.  From  this  one  may 
question  the  validity  of  the  normal  mode  method  for  evaporating  jets.  Another  possible  reason 
for  the  differences  can  be  tracked  to  their  initial  conditions  where  at  /=0  the  free  surface  is 
moving  but  there  is  no  internal  flow  field.  We  could  not  produce  such  initial  condition  in  our 
numerical  simulations. 


3.3  Breakup  time  and  volume 

Linear  theory  predicts  that  the  minimum  breakup  time  occurs  at  k=0.691,  which 
corresponds  to  the  wavenumber  with  the  maximum  growth  rate.  However,  this  is  only  true  for 
£0— >0,  otherwise  the  minimum  breakup  time  occurs  at  higher  values  of  k  as  explained  by 
Chaudhary  and  Redekopp  (1980).  In  the  simulations  presented  in  Fig.  7a,b,c,d  the  initial 
disturbance  amplitude  equals  0.05  therefore  most  of  our  curves  give  a  minimum  breakup  time  for 
&=0.8.  From  Fig.  7a  we  can  quantitatively  deduce  that  as  we  increase  the  evaporation  rate  we 
decrease  the  breakup  time,  which  has  been  demonstrated  by  other  graphs  presented  earlier.  This 
is  due  to  the  fact  that  the  addition  of  recoil  force  has  a  destabilizing  effect  for  low  speed  jets, 
which  we  are  considering  in  this  study.  Also  evaporation  decreases  the  mean  radius  of  our  jet 
with  time.  Figure  7b  clearly  shows  that  as  we  decrease  the  Reynolds  number  breakup  time 
increases,  which  is  due  to  the  additional  viscous  damping.  This  result  is  also  in  agreement  with 
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Figure  7:  Breakup  time  of  jet:  (a)  Re=200,zo=0.05,s=l0~3;  (b)  eo=0.05,p=3xl0"2,5=10"3;  (c) 
^=200,(3=3xl0"2,5=10"3;  (d)  /?e=200,£o=0.05,(5=3xl0"2. 


previous  observation  for  non-evaporating  jets  in  the  Reyleigh  regime  (Mashayek  and  Ashgriz 
1995).  Another  important  point  to  note  is  that  the  minimum  breakup  time  shifts  to  lower 
wavenumbers  as  Re  is  significantly  reduced.  This  is  in  accordance  to  linear  theory  which  states 
that  the  maximum  growth  rate  shifts  toward  the  lower  wavenumbers  for  more  viscous  liquid  jets. 
Even  though  our  cases  involve  evaporation,  Fig.  7c  still  re-emphasizes  Chaudhary's  results  by 
showing  that  as  e0  increases  from  zero  the  minimum  breakup  time  occurs  at  higher 
wavenumbers.    For  example,  for  eo=0.5  the  minimum  breakup  time  occurs  at  k=0.9  but  for 
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£o=0.001  this  happens  at  £=0.8.  Clearly,  from  Fig.  7c  we  notice  as  initial  disturbance  amplitude 
increases  breakup  time  decreases.  From  inspection  of  Fig.  7d  we  can  dissimulate  that  an 
increase  in  density  ratio  increases  breakup  time  due  to  the  diminishing  of  the  recoil  force. 

Breakup  volume,  the  instantaneous  combined  volume  of  both  the  main  drop  and  the 
satellite  drop  at  breakup,  is  primarily  dependent  on  jet  breakup  time.  If  it  takes  longer  for  jet 
breakup  to  occur  then  the  jet  is  allowed  more  time  to  evaporate.  Beginning  with  an  initial 
normalized  jet  volume  of  one,  Fig.  8  shows  the  final  volume  of  the  jet  while  varying  different 
parameters.  Analysis  of  the  results  from  each  graph  gives  similar  conclusions  to  that  described 
for  the  breakup  time. 


2 

CQ 


(a) 

1 


■  P=1x10 
♦  (3=2x10"' 
▲  0=3x10" 


(b) 

1 


MRe=W 
A  Re=20O 


o 

I 

pa 


(c) 

1 


0.2  0.4 


AE0=o.noi 

♦  E„=0-05 
■  e„=0.2 

•  c=0.5 


(d) 

1 


■ 

■ 

■ 

■ 

■ 

■ 

■ 

A 

A 

A 

A 

A 

A 

A 
♦ 

♦ 

♦ 

* 

♦ 

♦ 

♦ 

■  s=lxW 
As=lxl0~' 

+  s=lxl0~' 

Wavenumber 


Wavenumber 


Figure  8:  Breakup  volume  of  jet:  (a)  Re=200,Eo=0.05,s=10~3;  (b)  eo=0.05,P=3xl0"2,5=10"3;  (c) 
/?<?=200,p=3xl0-2,s=lfJ3;  (d)  /te=200,£o=0.05,p=3xlfJ2. 
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3.4  Drop  size  at  breakup 

The  control  of  satellite  drops  has  been  one  of  the  primary  driving  forces  behind  nonlinear 
jet  instability  studies.  For  non-evaporating  jets  the  total  volume  of  a  single  wavelength  of  fluid 
is  constant  throughout  the  breakup  process  and  equal  to  the  combination  of  the  main  drop 
volume  and  satellite  drop  volume  at  breakup.  However  for  evaporating  jets  the  volume  changes 
in  time  and  the  total  volume  of  the  jet  at  any  instant  is  equal  to  the  initial  jet  volume  minus  the 
total  amount  of  evaporation  up  to  that  point.  At  breakup,  the  volume  of  an  evaporating  jet  is 
equal  to  the  combined  volume  of  both  the  main  and  satellite  drops,  where  this  total  could  be 
significantly  less  than  the  initial  volume  depending  on  the  evaporation  rate  and  the  breakup  time. 
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Figure  9:  Variation  of  drop  size  at  breakup,  /te=200,eo=0.05,.s=10" 
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Figure  9  plots  the  main  and  satellite  drop  radii  versus  wavenumber  for  (3=102,2xlO"2,3xl02. 
From  this  figure  we  see  that,  for  small  evaporation  rate,  as  wavelength  increases  the  drop  size 
generally  increases,  which  is  similar  to  non-evaporating  jets.  However,  large  evaporation  rates 
and  small  wavenumbers  (i.e.  £=0.3  for  P=2xl0~  and  £=0.3,0.4  for  P=3xl0~2)  we  notice  a 
decrease  in  both  main  and  satellite  drop  size  due  to  the  increased  breakup  time,  which  allows 
more  vapor  to  leave  the  surface  and  reduce  the  overall  jet  volume.  Two  general  conclusions 
from  Fig.  9  are:  (i)  as  evaporation  rate  increases  drop  size  decreases  for  all  disturbance 
wavenumbers,  and  (ii)  comparing  the  percent  reduction  in  drop  size  evaporation  affects  the 
satellite  drop  more  than  the  main  drop  (e.g.  by  increasing  (3  by  a  factor  of  3  at  £=0.3  the  main 
drop  radius  reduces  38%  and  the  satellite  drop  radius  reduces  50%,  or  at  £=0.9  this  reduces  the 
main  drop  by  10%  and  the  satellite  drop  by  29%). 


3.5  Evaporation  flux  on  free  surface 

In  our  evaporation  model  we  approximate  the  unit  area  evaporation  rate  of  our  perturbed 
jet  surface  with  that  of  a  sphere  with  the  same  local  radius  of  curvature.  Curvature  of  the  jet 
surface  is  determined  by  the  height  function  using  (15).  To  be  in  agreement  with  (9)  the 
curvature  of  the  surface  was  monitored  throughout  simulations  to  prohibit  negative  values  in 
evaporation  cases.  Figure  10  shows  the  curvature  of  the  jet  versus  the  axial  coordinate  z  for  the 
simulations  presented  in  Fig.  2.  In  every  case  the  curvature  increases  in  the  neck  region  as  we 
approach  the  breakup  time.  This  translates  to  an  increase  in  the  evaporation  mass  flux  as  shown 
in  (26).  That  is,  the  evaporation  flux  at  the  neck  is  higher  than  that  at  the  swell,  leading  to 
further  pinching  of  the  neck.  Also  the  recoil  force  is  stronger  in  the  neck  region  than  that  at  the 
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Figure  10:  Variation  of  the  curvature  of  the  jet,  fo=200,  eo=0.5,  s=10"3:  (a)  k=0.1,  [5=0;  (b) 
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crest,  causing  the  liquid  to  be  squeezed  into  the  crests  and  thus  accelerating  the  growth  of  the 
disturbance  at  the  neck.     The  point  of  highest  curvature  in  each  set  of  data  represents  the 
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instantaneous  minimum  radius  of  curvature.  At  the  breakup  time  we  see  a  drastic  increase  in  the 
curvature  which  is  associated  with  the  breakup  point.  By  investigating  Fig.  10a,c,e  we  can  easily 
see  that  the  breakup  point  is  closer  to  the  neck  for  increasing  evaporation  rates. 


3.6  Pressure  and  velocity  distributions 

To  provide  further  insight  into  the  dynamics  of  the  jet  the  velocity  and  pressure 
distributions  within  the  liquid  jet  have  been  investigated.  Figure  1 1  shows  the  pressure 
distribution  and  velocity  vectors  at  three  different  non-dimensional  times  for  a  non-evaporating 
jet.  At  all  the  times  shown  we  can  easily  see  the  flow  of  the  fluid  from  the  neck  to  the  swell. 
With  increasing  surface  deformation  a  significant  increase  of  velocity  in  the  neck  region  takes 
place.  Initially  the  pressure  distribution  is  uniform  within  the  jet  but  as  the  velocity  field 
develops  we  can  clearly  identify  areas  of  low  and  high  pressure.  At  the  final  time  shown,  which 
is  near  the  breakup  time,  we  can  distinguish  a  high-pressure  region  at  the  top  of  the  swell  due  to 
the  conversion  of  velocity  to  pressure.  The  middle  ligament  exhibits  a  lower  pressure  where  the 
velocities  are  higher.  For  comparison  Fig.  12  gives  the  results  for  an  evaporating  jet  with  similar 
parameters  as  Fig.  11.  From  the  legend  we  see  that  the  evaporating  case  produces  a  higher 
overall  pressure,  which  is  due  to  the  contribution  of  both  capillary  and  recoil  forces.  We  can 
clearly  identify  the  recoil  force  in  both  the  middle  ligament  and  the  main  drop  of  our  jet  where 
we  see  a  higher  pressure  along  the  surface.  It  should  be  mentioned  that  in  Figs.  1 1  and  12  the 
velocity  vectors  have  been  scaled  in  each  of  the  3  non-dimensional  times  for  clarity. 

By  increasing  the  disturbance  wavelength  and  decreasing  the  initial  disturbance 
amplitude  our  code  has  been  able  to  predict  multiple  satellite  drops,  which  has  been  investigated 
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Figure  1 1 :  Pressure  distribution  and  velocity  vectors  for  a  non-evaporating  jet,  Re=200,  k=0J, 
eo=0.05,  p=0,  s=\0'3. 


experimentally  by  Vassallo  and  Ashgriz  (1991).     Figures   13  and   14  give  insight  into  the 
differences  seen  between  non-evaporating  and  evaporating  jets  under  these  conditions.  From  the 
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Figure  12:  Pressure  distribution  and  velocity  vectors  for  an  evaporating  jet,  Re=200,  k=0.1, 
8o=0.05,  (3=3xl0"2,5=103. 


non-evaporating  case,  Fig.  13,  we  notice  that  nonlinear  effects  have  drastically  influenced  the 
surface  shape  near  breakup.  Here  we  can  identify  a  main  drop  and,  possibly,  multiple  satellite 
drops  which  after  breakup  would  either  combine  and  oscillate  or  separate.  Investigating  Fig.  14 
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Figure  13:  Formation  of  multiple  satellite  drops  for  a  non-evaporating  jet,  Re=200,  k=0.3, 
eo=0.001,p=0,5=10"3. 


Figure  14:  Formation  of  multiple  satellite  drops  for  an  evaporating  jet,  Re=200,  k=0.3,  eo=0.001, 
(3=10-3,5=103. 


for  evaporating  case,  we  see  a  different  surface  shape  due  to  the  relocation  of  the  breakup  point 
close  to  the  initial  neck.  This  jet  creates  one  small  liquid  ligament  and  two  larger  bulged 
volumes  of  liquid,  which  would  continue  to  interact  after  breakup. 


4.  Conclusion 

Numerical  simulation  is  used  to  investigate  the  instability  of  evaporating  capillary  jets. 
The  evaporation  model  is  based  on  the  dMaw,  widely  used  for  spherical  drops.  To  implement 
this  model  for  perturbed  jets,  the  conventional  d2-law  is  modified  by  substituting  the  local  radius 
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of  the  surface  for  the  radius  of  a  sphere.  The  need  for  a  solution  for  the  ambient  gas  is 
eliminated  considering  small  density  ratios.  Numerical  simulations  are  carried  out  using  a 
Galerkin  finite  element  method  with  penalty  function.  The  motion  of  the  free  surface  is  tracked 
using  the  Height-Flux  Method  developed  by  Mashayek  and  Ashgriz  (1993).  With  this  technique, 
the  evaporation  at  the  surface  is  easily  handled  by  modifying  the  amount  of  fluid  within  each 
subvolume.  Preliminary  simulations  were  performed  to  establish  the  required  mesh  size  and 
time  increment  for  various  cases. 

In  non-dimensional  form,  the  formulation  involves  five  parameters:  the  amplitude  of  the 
initial  surface  disturbance,  the  wavenumber  of  the  initial  surface  disturbance,  the  Reynolds 
number,  the  non-dimensional  evaporation  rate,  and  the  density  ratio.  A  wide  range  of  variation 
for  each  of  these  parameters  is  considered  in  numerical  simulation.  The  results  were  discussed 
for  various  cases  leading  to  the  foremost  conclusion  that  evaporation  clearly  has  an  effect  on  jet 
instability,  mainly  increasing  the  growth  rate  of  the  disturbances  by  increasing  the  evaporation 
flux  at  the  neck.  Thus,  evaporation  causes  the  jet  to  break  up  sooner  and  reduces  the  size  of  both 
main  and  satellite  drops.  An  increase  in  internal  pressure  is  also  noticed  with  evaporating  cases 
due  to  the  contribution  of  both  capillary  and  recoil  forces.  Due  to  the  limited  number  of 
evaporating  jet  instability  studies  the  majority  of  our  conclusions  could  not  be  compared  to 
experimental  or  analytical  results. 
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